// Sample.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"

#include <stdexcept>
#include <iostream>
#include <map>
#include <ctime>
#include <cstdlib>
#include "string.h"
#include "blinalg.hpp"
#include "superelement.hpp"
#include "mesh.hpp"
#include "gauss.hpp"
#include "Intergration.hpp"
#include "functions.hpp"
#include "compute_matrix.h"
#define Max 50


int _tmain(int argc, _TCHAR* argv[])
{
	myindex N = 3;
	int max_refine = 2;
	std::cout << "Solve the 1D Elliptic problems: -(pu')'+qu=f by FEM" << std::endl;
	
	myindex nr_elements;
//	read and refine mesh
	nr_elements = readMeshFE("fem");
	for(int i=0;i<max_refine;i++)
		nr_elements *= 2;

	N = nr_elements - 1;
//	initialize superelements
	std::vector<FElement> elem(nr_elements);
	InitialMeshFE(elem,"fem");
	RefineMeshFE (elem, max_refine);

	std::cout << "Computing B..." << std::endl;
	FullMatrix<double> StiffnessMatrix(N,N);
	compute_matrix compute(N);
	compute.getStiffnessMatrix(StiffnessMatrix, elem);

	std::cout << "Computing F..." << std::endl;	
	//Tinh ve phai
	double d[2];
	//double *vephai = new double[N];
	//compute.getFVector(vephai, elem);
	EVector<double> vephai(N),ve(N);
	vephai.setZero();
	ve.setZero();

	for(int i=0;i<elem.size();i++)
	{
		Intergration(elem[i],d,10);
		if(i>0)
		{
			vephai[elem[i].getVertexIndex(0)-2] += d[0];
			ve[i-1] += d[0];
		}
		if(i<N-1)
		{
			vephai[elem[i].getVertexIndex(1)-2] += d[1];
			ve[i] += d[1];
		}
	}
	vephai.print();
	ve.print();
	return 0;
}

